Part 1では、データセット全体から「T細胞」という大きなグループを特定し、抽出しました。 しかし、一口にT細胞と言っても、その中には免疫を司令するヘルパーT細胞(CD4+)、攻撃を担うキラーT細胞(CD8+)、過剰な免疫を抑える制御性T細胞(Treg)など、多様な役割を持つ細胞が混ざっています。
全体解析で見つかった「高変動遺伝子(Variable Features)」は、T細胞とB細胞やマクロファージを分けるための遺伝子が中心です。 T細胞の内部にあるわずかな違いを浮き彫りにするためには、T細胞の中だけで遺伝子のばらつきを計算し直す必要があります。
前回保存したT細胞のサブセットデータを読み込みます。
# 保存したRDSファイルを読み込み
t_cells <- readRDS("./JIAsyno_CITEseq_T_share.rds")
# 読み込んだデータの確認(細胞数やAssayの種類を確認)
t_cells
## An object of class Seurat
## 30940 features across 11058 samples within 2 assays
## Active assay: RNA (30803 features, 2000 variable features)
## 3 layers present: counts, data, scale.data
## 1 other assay present: ADT
## 2 dimensional reductions calculated: pca, umap
T細胞集団に対して、改めて「特徴探し」から「分類」までを行います。
T細胞のサブタイプ(CD4, CD8, Tregなど)を区別するために重要な遺伝子を新たに選びます。
# ターゲットをRNAに設定
DefaultAssay(t_cells) <- "RNA"
# T細胞の中での高変動遺伝子を特定
t_cells <- FindVariableFeatures(t_cells, selection.method = "vst", nfeatures = 2000)
# 上位の変動遺伝子をプロット
top10_t <- head(VariableFeatures(t_cells), 10)
LabelPoints(plot = VariableFeaturePlot(t_cells), points = top10_t, repel = TRUE)
T細胞内の差異に基づいた新しい「軸(主成分)」を作成します。
t_cells <- ScaleData(t_cells)
t_cells <- RunPCA(t_cells, features = VariableFeatures(object = t_cells))
# 実際の主成分をプロットして確認
## ここでの色はT細胞以外の細胞を含めてクラスタリングしたときの情報です
DimPlot(t_cells, reduction = "pca") + NoLegend()
# どの主成分までが重要か、エルボープロットで確認
ElbowPlot(t_cells)
T細胞の中での「似た者同士」をグループ化します。
# PC1から10までを使用して計算
t_cells <- FindNeighbors(t_cells, dims = 1:10)
# resolution(解像度)を調整することで、クラスターの細かさを変えられます
# 今回は少し細かめに 0.8 に設定してみます
t_cells <- FindClusters(t_cells, resolution = 0.8)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11058
## Number of edges: 365998
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8853
## Number of communities: 18
## Elapsed time: 0 seconds
# UMAPの計算
t_cells <- RunUMAP(t_cells, dims = 1:10)
# 結果の描画
DimPlot(t_cells, reduction = "umap", label = TRUE)
各クラスターがどのようなT細胞なのか、代表的なマーカー遺伝子の発現を見てみましょう。
T細胞アノテーションのコツ まず、T細胞を大きく「CD4陽性」と「CD8陽性」に分け、その後にTregやMAIT、NKTなどの特殊なサブタイプを見ていくと良いでしょう。
なぜNK細胞のマーカーを見るのか? NK細胞はT細胞と共通の前駆細胞から分化し、かつ遺伝子発現状態が連続的につながっているため、T細胞集団に混入していることがあります。特にCD8陽性T細胞とNK細胞は一部マーカーが重複するため、注意が必要です。
# マーカー遺伝子の発現を地図(UMAP)上にプロット
features_to_plot <- c("CD3E", "CD4", "CD8A", "CCR7", "FOXP3", "KLRB1",
"TRAV1-2", "SLC4A10", "CXCL13",
"NCAM1", "FCGR3A",
"GZMB", "PRF1", "GZMK",
"TRDV2", "TRGV9",
"ZNF683", "KLRC2")
FeaturePlot(t_cells, features = features_to_plot, ncol = 4)
数値的な違いをより明確に確認します。
VlnPlot(t_cells, features = features_to_plot, pt.size = 0, ncol = 2)
既存のマーカーの発現以外に、データドリブンで各クラスターの特徴的な遺伝子を調べることも有効です。
t_markers = FindAllMarkers(t_cells, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
t_markers %>%
dplyr::group_by(cluster) %>%
dplyr::top_n(n = 8, wt = avg_log2FC) %>%
as.data.frame()
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 1.452093e-249 1.488029 0.807 0.480 4.472882e-245 0 NOSIP
## 2 1.361467e-230 1.544624 0.669 0.286 4.193727e-226 0 CCR7
## 3 5.819414e-202 1.796984 0.539 0.211 1.792554e-197 0 MAL
## 4 2.072956e-176 1.918416 0.408 0.130 6.385325e-172 0 FHIT
## 5 1.745770e-120 1.561202 0.470 0.225 5.377494e-116 0 MYC
## 6 1.785822e-113 1.781299 0.354 0.138 5.500868e-109 0 TRABD2A
## 7 6.781884e-71 1.359480 0.399 0.222 2.089024e-66 0 AC243960.1
## 8 6.444095e-69 1.673788 0.267 0.112 1.984975e-64 0 ACTN1
## 9 8.387867e-146 1.123077 0.924 0.694 2.583715e-141 1 LTB
## 10 5.609748e-132 1.623222 0.526 0.213 1.727971e-127 1 AQP3
## 11 4.280458e-108 1.348620 0.408 0.142 1.318509e-103 1 CD4
## 12 5.974422e-92 1.406259 0.459 0.202 1.840301e-87 1 TRADD
## 13 1.612873e-86 1.100976 0.490 0.224 4.968134e-82 1 TIMP1
## 14 5.828424e-56 1.286947 0.327 0.147 1.795329e-51 1 TNFRSF25
## 15 2.003784e-49 1.183526 0.367 0.186 6.172255e-45 1 LIMS1
## 16 6.814259e-38 1.180272 0.268 0.129 2.098996e-33 1 PASK
## 17 0.000000e+00 2.657202 0.841 0.243 0.000000e+00 2 CD8A
## 18 0.000000e+00 2.373755 0.784 0.217 0.000000e+00 2 CD8B
## 19 0.000000e+00 2.340439 0.783 0.243 0.000000e+00 2 GZMK
## 20 0.000000e+00 2.383701 0.658 0.157 0.000000e+00 2 GZMH
## 21 1.200088e-298 2.750063 0.556 0.116 3.696631e-294 2 LAG3
## 22 4.793321e-190 2.359342 0.415 0.094 1.476487e-185 2 HLA-DRB5
## 23 2.051888e-178 2.249384 0.437 0.107 6.320432e-174 2 HLA-DRA
## 24 3.221928e-160 2.186696 0.451 0.126 9.924506e-156 2 PTMS
## 25 0.000000e+00 5.237402 0.520 0.018 0.000000e+00 3 TRDV2
## 26 0.000000e+00 3.783232 0.504 0.061 0.000000e+00 3 TRGV9
## 27 4.069359e-254 2.453918 0.592 0.156 1.253485e-249 3 KLRG1
## 28 1.386016e-158 2.206775 0.307 0.057 4.269344e-154 3 ZBTB16
## 29 2.446471e-147 1.782404 0.414 0.108 7.535865e-143 3 PLEK
## 30 1.200782e-134 1.993738 0.415 0.119 3.698769e-130 3 MYBL1
## 31 5.954040e-122 1.817187 0.258 0.050 1.834023e-117 3 S1PR5
## 32 3.238299e-89 1.805637 0.291 0.083 9.974933e-85 3 C1orf21
## 33 0.000000e+00 2.485121 0.930 0.214 0.000000e+00 4 CD8B
## 34 0.000000e+00 3.223103 0.718 0.092 0.000000e+00 4 LINC02446
## 35 0.000000e+00 4.140657 0.294 0.016 0.000000e+00 4 NT5E
## 36 3.315425e-295 2.757895 0.537 0.087 1.021250e-290 4 AIF1
## 37 1.267340e-230 2.279822 0.511 0.105 3.903786e-226 4 ACTN1
## 38 1.840290e-218 2.011196 0.618 0.165 5.668644e-214 4 NELL2
## 39 4.573007e-187 2.799118 0.284 0.037 1.408623e-182 4 LRRN3
## 40 5.622938e-170 1.848587 0.529 0.141 1.732033e-165 4 TRABD2A
## 41 0.000000e+00 6.034507 0.807 0.032 0.000000e+00 5 FCGR3A
## 42 0.000000e+00 5.856208 0.834 0.062 0.000000e+00 5 SPON2
## 43 0.000000e+00 5.438919 0.814 0.048 0.000000e+00 5 FGFBP2
## 44 0.000000e+00 4.899510 0.627 0.031 0.000000e+00 5 KLRF1
## 45 0.000000e+00 5.562712 0.496 0.019 0.000000e+00 5 CX3CR1
## 46 0.000000e+00 5.246828 0.439 0.019 0.000000e+00 5 PRSS23
## 47 0.000000e+00 5.181054 0.383 0.021 0.000000e+00 5 AKR1C3
## 48 0.000000e+00 5.311373 0.345 0.022 0.000000e+00 5 MYOM2
## 49 1.877812e-277 3.307008 0.444 0.058 5.784226e-273 6 FXYD2
## 50 1.115310e-269 2.826102 0.609 0.122 3.435490e-265 6 ZNF683
## 51 5.236125e-229 3.638342 0.312 0.033 1.612884e-224 6 TRDV1
## 52 5.900723e-161 2.485535 0.485 0.117 1.817600e-156 6 KLRC2
## 53 1.769435e-105 1.983293 0.420 0.119 5.450391e-101 6 LINC02446
## 54 7.206443e-73 1.986776 0.335 0.105 2.219801e-68 6 RRAS2
## 55 2.375158e-67 1.976763 0.345 0.120 7.316199e-63 6 SPINT2
## 56 5.472146e-60 2.003276 0.273 0.082 1.685585e-55 6 IFNG-AS1
## 57 0.000000e+00 2.709259 0.854 0.197 0.000000e+00 7 FOSB
## 58 3.217484e-296 2.410720 0.989 0.405 9.910815e-292 7 FOS
## 59 4.161942e-257 2.377660 0.816 0.249 1.282003e-252 7 CREM
## 60 3.515260e-230 2.321351 0.951 0.529 1.082805e-225 7 TNFAIP3
## 61 1.071982e-180 2.893513 0.369 0.055 3.302026e-176 7 CDKN1A
## 62 9.548864e-177 2.367896 0.519 0.120 2.941337e-172 7 CSRNP1
## 63 1.406000e-149 2.337080 0.476 0.118 4.330902e-145 7 RGS2
## 64 3.060242e-127 2.767755 0.283 0.048 9.426463e-123 7 ATF3
## 65 0.000000e+00 2.684021 0.744 0.137 0.000000e+00 8 FHIT
## 66 2.186269e-228 2.215965 0.559 0.104 6.734364e-224 8 PLCL1
## 67 4.853353e-222 3.388579 0.289 0.026 1.494978e-217 8 AL589693.1
## 68 6.584765e-209 2.524472 0.454 0.073 2.028305e-204 8 LEF1-AS1
## 69 1.827189e-194 2.846443 0.292 0.031 5.628290e-190 8 ANKRD55
## 70 1.208567e-186 2.386368 0.389 0.059 3.722749e-182 8 AK5
## 71 1.335474e-166 2.306940 0.353 0.054 4.113662e-162 8 ADTRP
## 72 9.660521e-144 2.575178 0.274 0.038 2.975730e-139 8 PCSK5
## 73 0.000000e+00 2.997910 0.937 0.137 0.000000e+00 9 TYROBP
## 74 0.000000e+00 3.198042 0.849 0.231 0.000000e+00 9 GNLY
## 75 0.000000e+00 3.835183 0.651 0.068 0.000000e+00 9 KLRC1
## 76 5.142313e-305 4.400241 0.311 0.020 1.583987e-300 9 KIR2DL4
## 77 1.118948e-304 3.864970 0.363 0.029 3.446696e-300 9 NCAM1
## 78 1.389849e-280 3.581162 0.398 0.041 4.281152e-276 9 ATP8B4
## 79 3.353852e-208 3.626918 0.298 0.030 1.033087e-203 9 ITGAX
## 80 6.436316e-136 3.092236 0.266 0.035 1.982578e-131 9 GSN
## 81 0.000000e+00 5.776996 0.518 0.022 0.000000e+00 10 CXCL13
## 82 0.000000e+00 4.218815 0.476 0.022 0.000000e+00 10 PTPN13
## 83 0.000000e+00 4.057884 0.420 0.031 0.000000e+00 10 LINC00892
## 84 3.034151e-276 3.275340 0.420 0.041 9.346097e-272 10 TOX2
## 85 1.768715e-251 3.608642 0.314 0.023 5.448173e-247 10 ETV7
## 86 2.685888e-228 3.611146 0.334 0.031 8.273341e-224 10 GADD45G
## 87 1.667778e-210 3.682664 0.258 0.019 5.137255e-206 10 MYO7A
## 88 6.280697e-149 3.718157 0.256 0.028 1.934643e-144 10 NMB
## 89 0.000000e+00 6.387240 0.755 0.016 0.000000e+00 11 FOXP3
## 90 0.000000e+00 4.916334 0.476 0.033 0.000000e+00 11 IL2RA
## 91 7.505495e-299 3.383085 0.685 0.114 2.311918e-294 11 CTLA4
## 92 4.622075e-237 4.181195 0.262 0.016 1.423738e-232 11 LAYN
## 93 3.440854e-224 3.458952 0.427 0.052 1.059886e-219 11 RTKN2
## 94 7.695887e-211 3.254710 0.410 0.051 2.370564e-206 11 NCF4
## 95 5.061275e-148 3.386722 0.273 0.031 1.559025e-143 11 F5
## 96 4.305182e-106 3.526452 0.366 0.081 1.326125e-101 11 HPGD
## 97 0.000000e+00 4.584770 0.459 0.029 0.000000e+00 12 TNFSF9
## 98 8.334722e-216 3.577250 0.802 0.208 2.567344e-211 12 CCL4
## 99 5.797459e-137 3.439870 0.402 0.069 1.785791e-132 12 IFNG
## 100 3.498851e-136 2.738785 0.451 0.086 1.077751e-131 12 NR4A3
## 101 1.571180e-134 3.895646 0.371 0.058 4.839706e-130 12 CCL4L2
## 102 6.390201e-127 2.868786 0.368 0.059 1.968374e-122 12 HLA-DQA1
## 103 8.836986e-90 2.939809 0.337 0.070 2.722057e-85 12 CRTAM
## 104 4.691698e-82 2.716458 0.322 0.070 1.445184e-77 12 IER5L
## 105 0.000000e+00 3.868298 0.684 0.059 0.000000e+00 13 FXYD2
## 106 2.689449e-275 2.847425 0.740 0.117 8.284310e-271 13 KLRC2
## 107 3.002981e-210 2.890208 0.663 0.127 9.250082e-206 13 AREG
## 108 5.477543e-201 3.655236 0.372 0.038 1.687247e-196 13 TRDV1
## 109 5.307458e-180 4.035854 0.276 0.022 1.634856e-175 13 KLF4
## 110 3.444784e-156 4.061209 0.265 0.023 1.061097e-151 13 PLK2
## 111 1.143226e-130 3.053326 0.265 0.029 3.521478e-126 13 SPRY2
## 112 1.792414e-125 3.390112 0.288 0.038 5.521172e-121 13 GEM
## 113 8.577664e-74 3.948375 0.302 0.054 2.642178e-69 14 AC016831.7
## 114 1.953122e-60 3.324848 0.617 0.311 6.016203e-56 14 CHST11
## 115 1.307790e-49 3.327553 0.510 0.228 4.028386e-45 14 ATXN1
## 116 1.893749e-39 3.317182 0.328 0.104 5.833314e-35 14 ZNF407
## 117 1.340319e-30 3.300858 0.282 0.092 4.128584e-26 14 ATG7
## 118 1.924731e-28 3.134479 0.315 0.121 5.928750e-24 14 TRPS1
## 119 1.658359e-27 3.307763 0.273 0.094 5.108243e-23 14 ZSWIM6
## 120 3.755792e-27 3.112201 0.347 0.150 1.156897e-22 14 MGAT5
## 121 2.632881e-247 3.797971 0.862 0.134 8.110063e-243 15 XCL1
## 122 7.243149e-237 5.088743 0.792 0.130 2.231107e-232 15 AREG
## 123 1.934716e-219 3.793470 0.712 0.100 5.959507e-215 15 NR4A1
## 124 1.894780e-206 4.061271 0.762 0.129 5.836490e-202 15 XCL2
## 125 5.436837e-195 4.431189 0.323 0.020 1.674709e-190 15 KRT86
## 126 2.253200e-134 3.853521 0.769 0.217 6.940532e-130 15 CCL4
## 127 4.438854e-103 4.470166 0.338 0.044 1.367300e-98 15 CCL3
## 128 8.504690e-68 3.685516 0.369 0.073 2.619700e-63 15 CD83
## 129 0.000000e+00 8.780390 0.751 0.003 0.000000e+00 16 RRM2
## 130 0.000000e+00 11.126834 0.574 0.001 0.000000e+00 16 UBE2C
## 131 0.000000e+00 9.926383 0.541 0.001 0.000000e+00 16 BIRC5
## 132 0.000000e+00 8.864625 0.440 0.001 0.000000e+00 16 MYBL2
## 133 0.000000e+00 10.851496 0.402 0.000 0.000000e+00 16 SPC25
## 134 0.000000e+00 8.822303 0.383 0.002 0.000000e+00 16 CCNB2
## 135 0.000000e+00 9.491492 0.340 0.000 0.000000e+00 16 HJURP
## 136 0.000000e+00 10.465019 0.311 0.000 0.000000e+00 16 DLGAP5
## 137 1.532190e-66 4.129490 0.750 0.040 4.719604e-62 17 MYOM2
## 138 2.889103e-29 2.345311 0.833 0.107 8.899303e-25 17 SPON2
## 139 4.506285e-28 2.434113 0.708 0.077 1.388071e-23 17 FCGR3A
## 140 2.254802e-20 3.088388 0.250 0.015 6.945467e-16 17 NMUR1
## 141 1.970125e-16 2.384756 0.417 0.047 6.068576e-12 17 CX3CR1
## 142 2.809336e-12 2.491478 0.333 0.040 8.653599e-08 17 AL589693.1
## 143 4.336779e-09 2.401853 0.333 0.047 1.335858e-04 17 GSN
## 144 5.470198e-06 2.338720 0.250 0.048 1.684985e-01 17 CYB561D2.1
DoHeatmap(t_cells, features = t_markers %>%
dplyr::group_by(cluster) %>%
dplyr::top_n(n = 8, wt = avg_log2FC) %>%
dplyr::pull(gene) %>%
unique())
解析の結果をもとに、クラスター番号を生物学的な名前に書き換えます。
論文内のアノテーションとの一致を確認します。
サブクラス解析を行うことで、全体の解析では一つの「T細胞」という塊だったものが、実は多様な集団の集まりであることが分かりました。 自己免疫疾患(今回はJIA)の炎症部位において、「どのT細胞サブセットが特に活性化しているのか」、あるいは「治療によってどの集団が減るのか」を詳細に追跡することが、病態理解の鍵となります。
resolution を 0.2 や 1.5 に変えると、クラスターの数はどう変化しますか?sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 26.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Tokyo
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.3 magrittr_2.0.3 Seurat_5.2.1 SeuratObject_5.0.2
## [5] sp_2.1-2 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.15.0 jsonlite_1.8.8
## [4] spatstat.utils_3.1-2 ggbeeswarm_0.7.2 magick_2.8.2
## [7] farver_2.1.1 rmarkdown_2.25 vctrs_0.6.5
## [10] ROCR_1.0-11 spatstat.explore_3.3-4 htmltools_0.5.7
## [13] sass_0.4.8 sctransform_0.4.1 parallelly_1.36.0
## [16] KernSmooth_2.23-22 bslib_0.6.1 htmlwidgets_1.6.4
## [19] ica_1.0-3 plyr_1.8.9 plotly_4.10.3
## [22] zoo_1.8-12 cachem_1.0.8 igraph_1.6.0
## [25] mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
## [28] Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
## [31] fitdistrplus_1.1-11 future_1.33.1 shiny_1.8.0
## [34] digest_0.6.33 colorspace_2.1-0 tensor_1.5
## [37] RSpectra_0.16-1 irlba_2.3.5.1 labeling_0.4.3
## [40] progressr_0.14.0 fansi_1.0.6 spatstat.sparse_3.1-0
## [43] httr_1.4.7 polyclip_1.10-6 abind_1.4-5
## [46] compiler_4.3.2 withr_2.5.2 fastDummies_1.7.3
## [49] highr_0.10 MASS_7.3-60 tools_4.3.2
## [52] vipor_0.4.7 lmtest_0.9-40 beeswarm_0.4.0
## [55] httpuv_1.6.13 future.apply_1.11.1 goftest_1.2-3
## [58] glue_1.6.2 nlme_3.1-163 promises_1.2.1
## [61] grid_4.3.2 Rtsne_0.17 cluster_2.1.4
## [64] reshape2_1.4.4 generics_0.1.3 gtable_0.3.4
## [67] spatstat.data_3.1-4 tidyr_1.3.0 data.table_1.16.0
## [70] utf8_1.2.4 spatstat.geom_3.3-5 RcppAnnoy_0.0.21
## [73] ggrepel_0.9.4 RANN_2.6.1 pillar_1.9.0
## [76] stringr_1.5.1 limma_3.58.1 spam_2.10-0
## [79] RcppHNSW_0.5.0 later_1.3.2 splines_4.3.2
## [82] dplyr_1.1.4 lattice_0.21-9 survival_3.5-7
## [85] deldir_2.0-2 tidyselect_1.2.0 miniUI_0.1.1.1
## [88] pbapply_1.7-2 knitr_1.45 gridExtra_2.3
## [91] bookdown_0.37 scattermore_1.2 xfun_0.41
## [94] statmod_1.5.0 matrixStats_1.2.0 stringi_1.8.3
## [97] lazyeval_0.2.2 yaml_2.3.8 evaluate_0.23
## [100] codetools_0.2-19 tibble_3.2.1 BiocManager_1.30.22
## [103] cli_3.6.2 uwot_0.1.16 xtable_1.8-4
## [106] reticulate_1.35.0 munsell_0.5.0 jquerylib_0.1.4
## [109] Rcpp_1.0.11 globals_0.16.2 spatstat.random_3.3-2
## [112] png_0.1-8 ggrastr_1.0.2 spatstat.univar_3.1-2
## [115] parallel_4.3.2 ellipsis_0.3.2 ggplot2_3.4.4
## [118] presto_1.0.0 dotCall64_1.1-1 listenv_0.9.0
## [121] viridisLite_0.4.2 scales_1.3.0 ggridges_0.5.5
## [124] purrr_1.0.2 rlang_1.1.2 cowplot_1.1.2
# mega biteで表示
object.size(t_cells) / 1000 / 1000
## 654.6 bytes